#include "field.h"

void field_gradient(Field *      f_grad_x,
                    Field *      f_grad_y,
                    const Field *f,
                    const double inv_dx,
                    const double inv_dy)
{
    CHECK_FIELD(f_grad_x);
    CHECK_FIELD(f_grad_y);
    CHECK_FIELD(f);
    CHECK_FIELDS(f_grad_x,f_grad_y);
    CHECK_FIELDS(f,f_grad_x);
    {
        double       ** gx    = ARRAY_OF(f_grad_x);
        double       ** gy    = ARRAY_OF(f_grad_y);
        const double ** a     = CONST_ARRAY_OF(f);
        const double    fac_x = 1.0 * inv_dx;
        const double    fac_y = 1.0 * inv_dy;
        const int       nx    = f->nx;
        const int       ny    = f->ny;
        const int       nx1   = nx-1;
        const int       ny1   = ny-1;
        int i,j;
        
        for(j=0;j<ny;++j)
        {
            //const int jm = (j <= 0)   ? ny1 : j-1;
            const int jp = (j >= ny1) ? 0   : j+1;
            
            for(i=0;i<nx;++i)
            {
              //  const int im = (i <= 0)   ? nx1 : i-1;
                const int ip = (i >= nx1) ? 0   : i+1;
                
                gx[j][i] = fac_x * ( a[j][ip] - a[j][i] );
                gy[j][i] = fac_y * ( a[jp][i] - a[j][i] );
                
            }
            
        }
        
    }
}

